#!/usr/bin/env python
# -*- coding: utf-8  -*-
#
# This code is free software: you can
# redistribute it and/or modify it under the terms of the GNU General Public
# License as published by the Free Software Foundation, version 2.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc., 51
# Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Copyright Izaskun Mallona
# izaskun.mallona@gmail.com

"""Recombination sites finder. Based on PNAS | April 27, 2010 | vol. 107 | no. 17 | 7847–7852
www.pnas.org/cgi/doi/10.1073/pnas.0908774107"""

import re, os
from Bio import SeqIO

recSites = {
    'Myers':r'CC[ACTG]CC[ACTG]T[ACTG][ACTG]CC[ACTG]C',
    'M26CRE' : r'ATGACGT',
    'CCAAT' : r'CCAAT',
    'oligoC' : r'CCCCGCA',
    'ade64095' : r'GGTCT[AG]GAC',
    'motif86' : r'[AT]TCGGCCGA'
    }

compiledRecSites = {}

#process summary
recDataFrame = []

#header has 13 elements
header=[
    "Opuntia_id", 
    "Myers",
    "M26CRE",
    "CCAAT",
    "oligoC",
    "ade64095",
    "motif86"
        ]

#header shoud be the first row in dataFrame
recDataFrame.insert(0,header)
    
for i in range(len(recSites)):
    compiledRecSites[recSites.keys()[i]] = re.compile(recSites.values()[i], re.IGNORECASE)

def searchRecSites(query):
    Myers = 'F'
    M26CRE = 'F'
    CCAAT = 'F'
    oligoC = 'F'
    ade64095 = 'F'
    motif86 = 'F'
  
    for i in range(len(compiledRecSites)):
        if (re.search(compiledRecSites['Myers'], query) is not None) :
            Myers= 'T'
        if (re.search(compiledRecSites['M26CRE'], query) is not None) :
            M26CRE= 'T'   
        if (re.search(compiledRecSites['CCAAT'], query) is not None) :
            CCAAT= 'T'
        if (re.search(compiledRecSites['oligoC'], query) is not None) :
            oligoC= 'T'  
        if (re.search(compiledRecSites['ade64095'], query) is not None) :
            ade64095= 'T'
        if (re.search(compiledRecSites['motif86'], query) is not None) :
            motif86= 'T'  
        
        row=[
            record.id, 
            Myers,
            M26CRE,
            CCAAT,
            oligoC,
            ade64095,
            motif86          
        ]
        
        return row

#opening the handle
handlePath = os.environ['HOME'] +"/ourAssemblyAnnotation" +"/uniprot_trembl.fasta"
handle = open(handlePath, "rU")

for record in SeqIO.parse(handle, "fasta"):
    recDataFrame.append(searchRecSites(record.seq))
    
handle.close()
                  
#write an archive with the recombination search results
RECDATAFRAME = "/scratch1/" + os.environ['PBS_JOBID'] + "/" + "recombinationDataFrame.csv"
recombinationDataFrame=open(RECDATAFRAME, "w")
csvWriter=csv.writer(recombinationDataFrame)
csvWriter.writerows(recDataFrame)
recombinationDataFrame.close()
